# -*- coding: utf-8 -*-
from sympy.solvers.solvers import solve
from sympy import symbols, pprint, roots
from mymatrix import print_mat 
def solve_system_linear_diff_eq(m):
    k, x, y, z = symbols('k x y z')
    msize = len(m)
    print msize
    for i in xrange(0, msize):
        m[i][i] = m[i][i] - k
    eq = 0
    if msize == 3:
        nstep = 3
    elif msize == 2:
        nstep = 1
    for i in xrange(0, nstep):
        eq1 = 1
        eq2 = 1
        for j in xrange(0, msize):
            eq1 = eq1 * m[j][(j+i)%msize]
            eq2 = eq2 * m[j][(-(j+1)+i)%msize]
        eq = eq + eq1 - eq2
    eq.simplify()
    print u"characteristic polynomial is"
    pprint(eq.expand())
    print
    rs = roots(eq, k)
    print "\nRoots: (root: multiplicity)\n",
    pprint(rs)
    mr = []
    for str in m:
        s = str[:]
        mr.append(s)
    for r in rs:
        eqs = []
        m = []
        for str in mr:
            s = str[:]
            m.append(s)
        for i in xrange(0, msize):
            m[i][i] = m[i][i].subs(k, r)
        if msize == 3:
            for str in m:
                eqs.append(str[0]*x+str[1]*y+str[2]*z)
        if msize == 2:
            for str in m:
                eqs.append(str[0]*x+str[1]*y)            
        print u"\nFor k = \"", r, "\" equation is "
        pprint(solve(eqs, [x, y, z]))
        print 
        print_mat(m)
        print
